#'
#' Kriging a continous variable into any specified intervals.
#' 
#'@description autokriging via functions in the sgeostat package
#'
#'@param env an env object created by env
#'@param scale target scale of the kriging  
#'@param plotdim the width and height of the whole plot  
#'@param num.lags,maxdist.pair,maxdist.krige,type.variogram initial parameters used in calculate the empirical variogram, see in details 
#'@param iter,range,nugget,still initial parameters used to estimate the theoretical viarogram
#'
#'@details 
#'  The entire process first requires the point function to create a point object, 
#'  the pair function to create a pair objects, 
#'  the est.variogram to create a variogram object, 
#'  and fit.variogram to fit a model to the variogram. 
#'
#' Fitting variogram was tricky and required small enough number of lags to succeed, eg.
#' num.lags < maxdist/mindist, and mindist=20.  
#' 
#' Initial parameter estimates had to be good, else no fit curve was drawn and colinearity reported. 
#' 
#' Suggested maximum distance parameter estimates are: 
#' maxdist.pair = shorter plot dimension + 1 (ex. for a 400 x 200 m plot, maxdist.pair=201) 
#' maxdist.krige = 61 for 20 m intervals (the longer the distance the more time kriging will take, 
#' but it should not be less than 41). 
#'
#' @return 
#' it returns an env with expected scale
#' 
#'@note 
#' This function require sgeostat library. 
#' 
#'
#'@author 
#' Written by Kyle Harms, modified by S. Lao. & Guochun Shen
#' 
#'@export
autokriging=function(env,scale=5,plotdim=c(400,400),num.lags=10,maxdist.pair=201,type.variogram="exponential",
                  maxdist.krige=31,iter=30,range=100,nugget=2,sill=20,plot.it=FALSE)
{
  require(sgeostat)
  
  elevorig=data.frame(x=env$x,y=env$y,elev=env$z)
  if(any(plotdim%%scale!=0))
    stop("The plot can not be equally divided by the gride size defined by the scale")
  
  elevpt=point(elevorig)
  elevpair=pair(elevpt,num.lags=num.lags,maxdist=maxdist.pair)
  elev.variogram=est.variogram(elevpt,elevpair,"elev")
  
  elev.variogram.mod<-sgeostat::fit.variogram(model=type.variogram,elev.variogram,nugget=nugget,
                                              sill=sill,range=range,plot.it=plot.it,iterations=iter)
  
  
  # now to create columns for the 5m intervals, example, for a 1000x500m plot, x goes from 0 to 200 and y from 0 to 100. 
  
  krige.x = scale*rep(0:(plotdim[1]/scale),rep(1+plotdim[2]/scale,1+plotdim[1]/scale))
  krige.y = scale*rep(0:(plotdim[2]/scale),1+plotdim[1]/scale)
  krige.to = data.frame(cbind(krige.x,krige.y))
  colnames(krige.to)=c("x","y")
  krige.pts=point(krige.to)
  
  krigedist=sgeostat::krige(krige.pts,elevpt,"elev",elev.variogram.mod,
                            maxdist=maxdist.krige,extrap=T)
  
  kr=env(x=krigedist$x,y=krigedist$y,z=round(krigedist$zhat,3))

  return(kr)
}
